#include "Astronomy/Celestial/Sun.hpp"
#include "Astronomy/CFK5.hpp"
#include "Astronomy/Celestial/Earth.hpp"
#include "Astronomy/Constant.hpp"
#include "Astronomy/Nutation.hpp"
#include "Astronomy/Precession.hpp"
#include "Astronomy/Transformation.hpp"

double CelestialSun::GeometricEclipticLongitude(const double T, int CalcItemsCount) noexcept
{
    /*地心黄经 = 日心黄经 + 180度*/
    // 180.0 °= PI rad;;
    const double value = CelestialEarth::EclipticLongitude(T, CalcItemsCount);
    // return Transformation::MapTo0To2PIRange(value);
    return value;
}

double CelestialSun::GeometricEclipticLatitude(const double T, int CalcItemsCount) noexcept
{
    /*地心黄纬 = -日心黄纬*/
    const double value = -CelestialEarth::EclipticLatitude(T, CalcItemsCount);
    return value;
    // return Transformation::MapTo0To2PIRange(value);
}

double CelestialSun::ApparentEclipticLongitude(double T, int CalcItemsCount) noexcept
{
    // 太阳地心黄道;
    double longitude       = CelestialSun::GeometricEclipticLongitude(T, CalcItemsCount);
    // double latitude = CelestialSun::GeometricEclipticLatitude(T, CalcItemsCount);
    // double R = CelestialEarth::RadiusVector(T, CalcItemsCount);

    double resultLongitude = longitude;

#if 0
    //修正精度;
    const double SunLonCorrection = CFK5::CorrectionInLongitude(longitude, latitude, T);

    resultLongitude = longitude + SunLonCorrection;
    
    // 修正天体章动;
    auto AdjustSunEclipticAngle = CNutation::HighPrecision(T, 0);

    resultLongitude += AdjustSunEclipticAngle.longitude;

    /*太阳地心黄经光行差修正*/
    //resultLongitude -= (0.005775518 * R * Transformation::DMSToDegrees(0, 0, VariationEclipticLongitudeHigh(T)));
    resultLongitude -= (0.005775518 * R *  VariationEclipticLongitudeHigh(T));

#else
    // 修正天体章动;
    resultLongitude += CNutation::MediumPrecisionLongitude(T);

    resultLongitude += VariationEclipticLongitude(T);

    /*[暂时不启用]太阳地心黄经光行差修正项是: -20".4898/R*/
    // resultLongitude -= Transformation::DMSToDegrees(0, 0, 20.4898 / R);
#endif

    return resultLongitude + Constants::PI;
}

double CelestialSun::ApparentEclipticLatitude(double T, int CalcItemsCount) noexcept
{
    // 太阳地心黄道;
    double longitude = CelestialSun::GeometricEclipticLongitude(T, CalcItemsCount);
    double latitude  = CelestialSun::GeometricEclipticLatitude(T, CalcItemsCount);
    latitude += CFK5::CorrectionInLatitude(longitude, T);
    return latitude;
}

Coordinate3D CelestialSun::ApparentEclipticLocation(const double T, bool bHighPrecision) noexcept
{
    // 太阳地心黄道;
    double longitude              = CelestialSun::GeometricEclipticLongitude(T, bHighPrecision); // 2.2608225268479778
    double latitude               = CelestialSun::GeometricEclipticLatitude(T, bHighPrecision);  // 1.5907354109501171e-06
    double R                      = CelestialEarth::RadiusVector(T, bHighPrecision);             // 1.0149471547041036

    // 修正精度;
    const double SunLonCorrection = CFK5::CorrectionInLongitude(longitude, latitude, T);
    const double SunLatCorrection = CFK5::CorrectionInLatitude(longitude, T);

    //;
    Coordinate3D location{longitude + SunLonCorrection, latitude + SunLatCorrection, R}; // 0.0000015907186312685412
    if (bHighPrecision)
    {
        // 修正天体章动;
        auto AdjustSunEclipticAngle = CNutation::HighPrecision(T, 0);

        location.X += AdjustSunEclipticAngle.longitude;

        /*太阳地心黄经光行差修正*/
        location.X -= (0.005775518 * R * Transformation::DMSToDegrees(0, 0, VariationEclipticLongitudeHigh(T)));
    }
    else
    {
        // 修正天体章动;
        location.X += CNutation::MediumPrecisionLongitude(T * 10.0);

        location.X += VariationEclipticLongitude(T);

        location.X = Transformation::MapTo0To2PIRange(location.X);

        /*[暂时不启用] 太阳地心黄经光行差修正项是: -20".4898/R*/
        // location.X -= Transformation::DMSToDegrees(0, 0, 20.4898/R);x = 2.26078530922819
    }
    return location;
}

double CelestialSun::VariationEclipticLongitudeHigh(const double T) noexcept
{
    const double Tsquared{T * T};
    const double Tcubed{Tsquared * T};

    const double deltaLambda{3548.193 + (118.568 * sin(Transformation::DegreesToRadians(87.5287 + (359993.7286 * T)))) +
                             (2.476 * sin(Transformation::DegreesToRadians(85.0561 + (719987.4571 * T)))) +
                             (1.376 * sin(Transformation::DegreesToRadians(27.8502 + (4452671.1152 * T)))) +
                             (0.119 * sin(Transformation::DegreesToRadians(73.1375 + (450368.8564 * T)))) +
                             (0.114 * sin(Transformation::DegreesToRadians(337.2264 + (329644.6718 * T)))) +
                             (0.086 * sin(Transformation::DegreesToRadians(222.5400 + (659289.3436 * T)))) +
                             (0.078 * sin(Transformation::DegreesToRadians(162.8136 + (9224659.7915 * T)))) +
                             (0.054 * sin(Transformation::DegreesToRadians(82.5823 + (1079981.1857 * T)))) +
                             (0.052 * sin(Transformation::DegreesToRadians(171.5189 + (225184.4282 * T)))) +
                             (0.034 * sin(Transformation::DegreesToRadians(30.3214 + (4092677.3866 * T)))) +
                             (0.033 * sin(Transformation::DegreesToRadians(119.8105 + (337181.4711 * T)))) +
                             (0.023 * sin(Transformation::DegreesToRadians(247.5418 + (299295.6151 * T)))) +
                             (0.023 * sin(Transformation::DegreesToRadians(325.1526 + (315559.5560 * T)))) +
                             (0.021 * sin(Transformation::DegreesToRadians(155.1241 + (675553.2846 * T)))) +
                             (7.311 * T * sin(Transformation::DegreesToRadians(333.4515 + (359993.7286 * T)))) +
                             (0.305 * T * sin(Transformation::DegreesToRadians(330.9814 + (719987.4571 * T)))) +
                             (0.010 * T * sin(Transformation::DegreesToRadians(328.5170 + (1079981.1857 * T)))) +
                             (0.309 * Tsquared * sin(Transformation::DegreesToRadians(241.4518 + (359993.7286 * T)))) +
                             (0.021 * Tsquared * sin(Transformation::DegreesToRadians(205.0482 + (719987.4571 * T)))) +
                             (0.004 * Tsquared * sin(Transformation::DegreesToRadians(297.8610 + (4452671.1152 * T)))) +
                             (0.010 * Tcubed * sin(Transformation::DegreesToRadians(154.7066 + (359993.7286 * T))))};
    return deltaLambda;
}

double CelestialSun::VariationEclipticLongitude(double T) noexcept
{
    const double Tsquared{T * T};

    double v = -0.043126 + 628.301955 * T - 0.000002732 * Tsquared; // 平近点角;
    double e = 0.016708634 - 0.000042037 * T - 0.0000001267 * Tsquared;

    return Transformation::ArcsecToRadians(-20.49552 * (1 + e * cos(v))); // 黄经光行差
}

double CelestialSun::RiseSet(double jd, double L, double fa, double sj)
{
    jd = floor(jd + 0.5) - L / Constants::PI2;
    for (int i = 0; i < 2; i++)
    {
        double T = jd / 36525.0;
        double E = Transformation::ArcsecToRadians(84381.1060 - 46.836769 * T);
        double J = (48950621.66 + 6283319653.318 * T + 53 * T * T - 994.0 + 334166.0 * cos(4.669257 + 628.307585 * T) + 3489.0 * cos(4.6261 + 1256.61517 * T) +
                    2060.6 * cos(2.67823 + 628.307585 * T) * T) /
                   10000000.0;
        // 太阳黄经以及它的正余弦值;
        double sinJ = sin(J), cosJ = cos(J);
        // 恒星时(子午圈位置);
        double gst =
            (0.7790572732640 + 1.00273781191135448 * jd) * Constants::PI2 + Transformation::ArcsecToRadians(0.014506 + 4612.15739966 * T + 1.39667721 * T * T);
        // 太阳赤经;
        double A     = atan2(sinJ * cos(E), cosJ);
        // 太阳赤纬;
        double D     = asin(sin(E) * sinJ);
        // 太阳在地平线上的cos(时角)计算;
        double cosH0 = Transformation::ArcsecToRadians(sin(-50.0 * 60.0) - sin(fa) * sin(D)) / (cos(fa) * cos(D));
        if (abs(cosH0) >= 1)
            return 0;
        //(升降时角-太阳时角)/太阳速度;
        jd += Transformation::MapToNegPIToPIRange(sj * acos(cosH0) - (gst + L - A)) / 6.28;
    }
    return jd;
}

double CelestialSun::GetEclipticLongitudeTime(double W) noexcept
{
    double t, v = 628.3319653318;
    t = (W - 1.75347 - Constants::PI) / v;
    v = CelestialEarth::Velocity(t); // v的精度0.03%，详见原文;
    t += (W - ApparentEclipticLongitude(t, 10)) / v;
    v = CelestialEarth::Velocity(t); // 再算一次v有助于提高精度,不算也可以;
    t += (W - ApparentEclipticLongitude(t, -1)) / v;
    return t;
}

double CelestialSun::GetEclipticLongitudeTimeLow(double W) noexcept
{
    double t, v = 628.3319653318;
    t = (W - 1.75347 - Constants::PI) / v;
    t -= (0.000005297 * t * t + 0.0334166 * cos(4.669257 + 628.307585 * t) + 0.0002061 * cos(2.67823 + 628.307585 * t) * t) / v;
    t += (W - CelestialEarth::EclipticLongitude(t, 8) - Constants::PI + Transformation::ArcsecToRadians((20.5 + 17.2 * sin(2.1824 - 33.75705 * t)))) / v;
    return t;
}
double CelestialSun::CalcTimeDifferenceLow(double t) noexcept
{
    double L = (1753470142.0 + 628331965331.8 * t + 5296.74 * t * t) / 1000000000.0 + Constants::PI;
    double E = Transformation::ArcsecToRadians(84381.4088 - 46.836051 * t);
    // 地球坐标;
    Coordinate3D z{CelestialEarth::EclipticLongitude(t, 5) + Constants::PI, 0, 0};

    // 太阳地心赤道坐标;
    z.EclipticToEquatorial(E);
    L = Transformation::MapToNegPIToPIRange(L - z.X);
    return L / Constants::PI2;
}
double CelestialSun::CalcTimeDifferenceHigh(double t) noexcept
{
    const double t2 = t * t;
    const double t3 = t2 * t;
    const double t4 = t3 * t;
    const double t5 = t4 * t;
    double       L  = (1753470142.0 + 628331965331.8 * t + 5296.74 * t2 + 0.432 * t3 - 0.1124 * t4 - 0.00009 * t5) / 1000000000.0 + Constants::PI -
               Transformation::ArcsecToRadians(20.5);
    // 黄经章;
    const double dL = Transformation::ArcsecToRadians(-17.2 * sin(2.1824 - 33.75705 * t)); // 黄经章
    // 交角章;
    const double dE = Transformation::ArcsecToRadians(9.2 * cos(2.1824 - 33.75705 * t)); // 交角章
    // 真黄赤交角;
    const double E  = Precession::EclipticObliquity(t) + dE;
    // 地球坐标;

    Coordinate3D z{CelestialEarth::EclipticLongitude(t, 50) + Constants::PI + VariationEclipticLongitude(t) + dL,
                   -(2796 * cos(3.1987 + 8433.46616 * t) + 1016 * cos(5.4225 + 550.75532 * t) + 804 * cos(3.88 + 522.3694 * t)) / 1000000000.0,
                   0};
    // 太阳地心赤道坐标;
    z.EclipticToEquatorial(E);
    z.X -= dL * cos(E);

    L = Transformation::MapToNegPIToPIRange(L - z.X);
    return L / Constants::PI2;
}